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As galaxy surveys begin to measure the imprint of baryonic acoustic oscillations (BAO) on large- 
scale structure at the sub-percent level, reconstruction techniques that reduce the contamination 
from nonlinear clustering become increasingly important. Inverting the nonlinear continuity equa¬ 
tion, we propose an Eulerian growth-shift reconstruction algorithm that does not require the dis¬ 
placement of any objects, which is needed for the standard Lagrangian BAO reconstruction al¬ 
gorithm. In real-space DM-only simulations the algorithm yields 95% of the BAO signal-to-noise 
obtained from standard reconstruction. The reconstructed power spectrum is obtained by adding 
specific simple 3- and 4-point statistics to the pre-reconstruction power spectrum, making it very 
transparent how additional BAO information from higher-point statistics is included in the power 
spectrum through the reconstruction process. Analytical models of the reconstructed density for 
the two algorithms agree at second order. Based on similar modeling efforts, we introduce four 
additional reconstruction algorithms and discuss their performance. 


I. INTRODUCTION 

Measuring the imprint of baryonic acoustic oscillations (BAO) on large-scale structure is becoming increasingly 
important for tightening constraints on the ACDM model parameters and probing the nature of dark energy and the 
geometry of the universe; see e.g. [I] for a recent review. The BAO lead to a bump in the 2-point halo correlation 
function at separation r ~ 150 Mpc, and to wiggles in the Fourier space power spectrumQ Measuring the BAO scale 
from this provides a standard ruler for a physical scale at low redshift which is very powerful in breaking parameter 
degeneracies that would be present when considering the high-redshift Cosmic Microwave Background (CMB) alone 
(e.g. [Ml). Measuring the distinct large-scale BAO bump in the correlation function or the corresponding wiggles in 
the power spectrum is also very robust against potential astrophysical and observational systematics. 

Nonlinear gravitational clustering smears the sharp linear BAO bump in the correlation function, mostly due to 
large-scale bulk flows and cluster formation [3]. In a procedure known as BAO reconstruction these large-scale shifts 
are reversed to (partially) recover the linear BAO signal j8]. Slightly more generally, the goal of large scale structure 
reconstruction is to restore information about the linear modes that were corrupted by nonlinear structure formation. 
In practice, BAO reconstruction works extraordinarily well and reduces the uncertainty on the BAO scale significantly, 
e.g. by a factor of more than 1.5 for BOSS DR11 [S] (see also nsHUD- Reconstruction also reduces a small bias of 
the BAO scale that is caused by nonlinear clustering. This has been studied in detail analytically and in simulations 

C5H23|. 

Although the standard Lagrangian reconstruction algorithm works very well, there are several motivations to explore 
new reconstruction algorithms, e.g. 

• to provide more insight where additional BAO information comes from, 

• derive algorithms more rigorously, 

• improve the methodology by requiring fewer model assumptions for data processing, 

• increase BAO information further. 

Motivated by this, we propose a number of new reconstruction algorithms. Inspired by the Lagrangian and Eulerian 
perspectives of fluid dynamics, where either individual particle positions or fields like the mass density are regarded 
as the fundamental dynamic variable, we separate these reconstruction algorithms into two broad, disjoint categories: 


1 This BAO imprint can be understood from a simple physical picture by considering perturbations to an initially point-like overdensity at 
the origin [2]. Dark matter just gravitationally clusters around the origin. Baryons and photons are tightly coupled before recombination 
(2 > 1100) and undergo acoustic oscillations, because infall due to gravitational clustering of baryons is counter-acted by photon pressure. 
As the universe cools, baryons and photons decouple, so that photons free-stream and baryons cluster in a spherical shell around the 
origin whose radius is given by the distance the sound wave traveled between the initial time and the time of decoupling; the BAO scale. 
The dark matter in the origin and the baryons in the shell around the origin subsequently cluster gravitationally and ultimately form 
halos which are separated by the BAO scale. For realistic initial overdensities, this effect is still present statistically. 
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Lagrangian reconstruction algorithms, which involve displacing individual objects at some point in the recon¬ 
struction process, and Eulerian reconstruction algorithms, which do not involve displacing individual objects at 
any point in the reconstruction process, e.g. by operating only on density fields. In this notion ‘mixed’ algorithms 
like the standard BAO reconstruction, which moves individual objects by a displacement field which is calculated 
from the linearized continuity equation for the mass density field, are still categorized as Lagrangian because they 
move individual objects at some step of the algorithm. A very simple (though not particularly useful) example for 
an Eulerian reconstruction algorithm is to suppress nonlinear peaks of the mass density by subtracting the squared 
density from the original density, i.e. <5 rec (x) = <5(x) — <5 2 (x), which does not involve displacing any objects. 

Although we cannot address all motivational bullet points above at once, we use them as guidelines for the Eulerian 
reconstructions introduced in this paper. Addressing the first point, we find that the power spectrum after our Eulerian 
reconstructions can be expressed as a sum of the unreconstructed power spectrum and specific simple 3- and 4-point 
statistics of the unreconstructed density. This makes it very transparent and explicit how the additional information on 
the BAO scale due to reconstruction comes from specific higher order N-point statistics of the unreconstructed density. 
Lagrangian reconstructions must also use higher order N-point information implicitly, because phase information of 
the density is used so that a reconstructed power spectrum realization cannot be obtained by any processing of a 
measured unreconstructed power spectrum alone. However, we are not aware of any way to write this in the same 
transparent and explicit form as for Eulerian reconstructions. 

Second, our Eulerian growth-shift reconstruction algorithm can be derived from the nonperturbative continuity 
equation that follows from mass conservation. It is obtained by expressing the time derivative of the nonlinear mass 
density in terms of the nonlinear mass density and the linearized velocity density. This is fully nonperturbative in the 
mass density, which is important because reconstruction is most interesting in the nonlinear regime where perturbation 
theory breaks down. We connect this Eulerian algorithm to the standard Lagrangian reconstruction algorithm by 
showing that these two algorithms are equivalent in second order perturbation theory. This theoretical connection 
can be regarded as a new argument for the robustness and success of the standard Lagrangian BAO reconstruction in 
the nonlinear regime; it automatically includes 3- and 4-point information in a specific form imposed by the nonlinear 
continuity equation. Note however that this point of view is limited by the fact that the nonlinear smoothing of 
the BAO is largely caused by third order and higher order perturbations, while we only show the equivalence of the 
algorithms up to second order. Still, small but important BAO shifts are caused by second order perturbations, and 
these shifts are reversed in a very similar fashion for Eulerian and Lagrangian reconstructions due to their equivalence 
at second order. 

Third, the methodology of Lagrangian reconstructions may be criticized for intermixing observed data with model 
assumptions that should actually be tested by the data. To see this, note that for Lagrangian reconstructions, the 
observed data is processed by changing observed galaxy positions according to a displacement field that is computed 
from the same observed galaxy positions under certain model assumptions (including assumptions on halo bias, 
redshift space distortions, or general relativity). Thus, the data and assumed fiducial model are mixed together 
before performing the likelihood analysis that compares the model-dependently-processed data against other models. 
Additionally, if the parameters for the displacement field depend on external datasets (e.g. the CMB to obtain the 
logarithmic growth rate through / ~ which holds in general relativity), then the reconstructed density is strictly 
speaking not independent from this external dataset any more, which may raise concerns about combining constraints 
on cosmological parameters from these datasets without potentially double-counting information. 

In practice, all of these concerns can of course be modeled or simulated by displacing objects with fiducial models 
that differ from the true simulated catalogs (e.g. different values of bias and logarithmic growth rate). The standard 
reconstruction algorithm seems to be rather robust to such changes. For example, reference m found that using a 
20%-high linear halo bias for the displacement field shifts the BAO scale by only ler or 2% for their experimental 
specifications. However, such shifts may be larger for other survey geometries and need to be checked for every new 
survey geometry m- Such tests should thus be an integral part of every reconstruction analysis, but running them for 
every survey geometry and BAO analysis renders the use of reconstruction somewhat cumbersome and computationally 
expensive. The Eulerian reconstructions proposed here can ameliorate these complications, because they only depend 
on the density field of the untransformed observed galaxy positions, and model-independent transformations of this 
field (like squaring it). Fiducial bias parameters simply enter the final reconstructed power spectrum by rescaling 
model-independently measured cross-spectra, so varying these parameters is computationally trivial. 

Finally, some of our Eulerian reconstruction algorithms could in principle increase BAO information further, because 
highly nonlinear scales are shifted differently than in standard reconstruction. In 2LPT, all reconstruction algorithms 
considered in this paper are equally well motivated, so it is not clear a priori clear which one recovers most BAO 
information. Comparing the performance of the algorithms in simulations is one of the main goals of this paper. 

The paper is organized as follows. We start by deriving two new Eulerian reconstruction schemes that work at 
the field level: First, the nonlinear continuity equation is used to derive the Eulerian growth-shift reconstruction 
algorithm. Second, we construct an Eulerian F -2 reconstruction algorithm that subtracts the second order F 2 part 
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of the density. We then put this into context by modeling the standard Lagrangian BAO reconstruction and some 
possible alternative Lagrangian reconstruction methods with 2LPT. Section |V| provides a high-level overview of the 
six reconstruction algorithms introduced in the previous sections. We then proceed with simulations and numerical 
results, discussing the origin of additional BAO information and the performance of all algorithms. Finally, we 
summarize and conclude. Appendices present Lagrangian modeling calculations, comparisons of reconstructions with 
2D slices, a more formal motivation of the Eulerian F 2 reconstruction using a Newton-Raphson iteration, and details 
of the numerical analysis. 


Notation and conventions 


We use the following Fourier conventior0 

/( k) = J d 3 x e- lk x /(x), /(x) 

i.e the Fourier transform of the gradient of a function is 


d 3 k 
(27r) 3 




7(k), 


[V/](k) = ik/(k). 


(1) 

( 2 ) 


Throughout the paper, q stands for Lagrangian configuration-space coordinates and x for Eulerian configuration- 
space coordinates, while k and k, denote Fourier-space wavevectors. Coordinates and wavenumbers are comoving. 
To shorten expressions we write 



d 3 k! d 3 k 2 
(27t) 3 (27t) 3 


(2^) 3 Mk 


ki 


k2)- 


(3) 


Throughout the paper, second order Lagrangian Perturbation Theory or 2LPT refers to expanding all terms in LPT 
consistently in the linear density S 0 up to 0(6q), i.e. we do not include corrections beyond second order that would 
follow e.g. by using the cumulant expansion theorem (which would be an interesting extension). 


II. EULERIAN GROWTH-SHIFT RECONSTRUCTION FROM NONLINEAR CONTINUITY 

EQUATION 

A. Dark Matter 

We propose to reverse the nonlinear time evolution of structures by starting with the observable nonlinear density 
S(x,ij) today and going back in time by some time interval Ar/ using a Taylor expansion in time, 

£( x , V — A 77 ) w <5(x, rj) - Ar] d v 6(x., rj) = S ie c (x), (4) 

where 77 is conformal time, 77 = f“ a f jq”') • The time derivative of the (nonlinear) density can be obtained from 
the nonperturbative continuity equation that follows from mass conservation. It states that any change of the mass 
density in a particular volume element with time must be caused by an equivalent in- or out-flow of mass (described 
by the divergence of the mass-weighted velocity or momentum), 

d v 5 + V • [(1 + £)v] = 0. (5) 

Here, v is the (nonlinear) peculiar velocity perturbation which we assume to be curl-free with 9 = V • v, i.e. v(k) = 
— zkd(k)/fc 2 . Writing out the spatial derivative gives 

<V(x, rj) = -V • v(x, 77 ) - v(x, rj) ■ V<5(x, rj) - S(x, r/)V • v(x, rj). (6) 

If we know the density and velocity fields then we can directly use this in Eq. Q to reconstruct the density at an 
earlier time. 


2 The sign in the exponent is the same as in m but opposite from |24| . 
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In practice, the locations of objects are easier to observe than their peculiar velocities, so that the mass density 
is typically much better known than the velocity field. At linear order we can infer the velocity field from the mass 
density, 

v o(k, 77) = fH^SoiK rj). (V 

This relation is obtained from the perturbative expansion 

S(x,r/) = D(t))6o(x) + 0(D 2 6$), (8) 

where time and spatial dependence are described by the linear growth function D(rj) and the linear perturbation 
<$o(x), respectively. The linearized continuity equation then becomes 

V • v 0 (x, 77 ) = -8 0 (-x.)d ri D(r)) = -fH6 0 (x, 77 ), (9) 

which implies Eq. ([7j) . Here, / = din D/d In a is the logarithmic growth rate and R = aH is the comoving Hubble 
parameter. The linear density $0 in Eq. ([T]) is unknown in observations, but it can be approximated by the observed 
nonlinear density if small scales corrupted by nonlinearities are smoothed out: 

5 0 (k)^S R (k) = W R (k)8(k), ( 10 ) 

where W R is a smoothing kernel that is unity on large scales k < R~ 1 and vanishes on small scales k > R-. 
The smoothing scale R should be chosen close to the scale where linear perturbation theory breaks down. The 
corresponding approximate velocity is 


v(x, 77 ) « —fRs(k, 77 ), 


( 11 ) 


where we defined 


s (k, 77 ) = -^W R (k)8(k, 77 ), 

which coincides with the negative Zeldovich displacement. The velocity divergence is simply 

V • v(x, 77) « - fUS R {x , 77 ). 

The density time derivative § then becomes 

d v 6(x,ri) = fH S R (x, 77 ) + s(x, 77 ) ■ V5(x, 77 ) + (5(x, r/)5 R (x, rj) , 
so that the Euler reconstructed field of Eq. 0 is given by 

5 rec (x) = 5(x, 77 ) - A 77 fR [s(x, 77 ) • V<5(x, 77 ) + S(x, r/)S R (x, 77 ) . 


( 12 ) 

(13) 

(14) 

( !5) 

We dropped the first term in the square brackets of Eq. ( |l4| ) because it describes linear time evolution but we only want 
to reverse nonlinear time evolution when reconstructing the linear density. The go-back time A 77 is a free parameter, 


A ?7 is 

but from now on we choose A 77 = (fR)- 1 so that the coefficient of the square brackets in Eq. ( |15| ) becomes —1. We 
call the resulting algorithm Eulerian growth-shift reconstruction, or short ‘EGS’, 

<5egs( x ) = <5(x, 77) - s(x, 77 ) • V<5(x, 77 ) - (5(x, r])S R {x, 77 ), 


(16) 


shift 


growth 


because it subtracts nonlinear growth and shift from the nonlinear density at the field level. As desired, the right 
hand side of Eq. (16) depends only on the mass density S and not the velocity so it can be easily obtained from 
observations. 

The shift and growth term in Eq. (16) are quadratic in the nonlinear unreconstructed density. The 2-point correlation 
function or power spectrum of the reconstructed density thus involves specific 2-, 3- and 4-point statistics of the 
unreconstructed density, 


-^EGS’ 5 BGs(k) = -^(k) ~^Ps 2 ,8(k) - 2P s .^s,s(k) + P S 2 S 2(k) + P s . vy s .v<5 (k) + 2P s . VSS 2(k), 


(17) 


2pt 


3pt 


4pt 
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where we sloppily wrote 1 5 2 for S R (x)S(x). This makes it very transparent how reconstruction adds information to the 
power spectrum by combining information from 2-, 3- and 4-point stat istics of the unreconstructed density. 

In fact, the 3-point cross-spectra of S 2 or s • V<5 with S in Eq. ( fl7| ) do not just probe the 3-point function, but 
they represent nearly-optimal maximum-likelihood estimators for the amplitudes of the contributions to the tree-level 
3-point function due to nonlinear gravitational growth and shift |24j . Roughly, they measure the projection of the 
observed 3-point function on tree-level expectations, and are thus expected to capture a significant fraction of the full 
3-point information. 

Note that the standard Lagrangian reconstruction method also employs the continuity equation. However, it is 
only used there in its fully linearized form, where both mass density and velocity are linearized, to obtain the first 
order velocity or displacement So = — ik/k 2 So. We also use this to approximate the velocity. In contrast to previous 
literature, we then plug this linearized velocity back into the nonlinear continuity equation without linearizing the 
mass density. The resulting density time derivative (141 is nonperturbatively correct in the unsmoothed nonlinear 
density (appearing on the left hand side and in the last two summands on the right hand side of Eq. @). It 
only uses perturbation theory to approximate the factors in Eq. ([6]) that involve the velocity, leading to s and 8 R 
in Eqs. (|14[) and (16). Thus, the crucial difference to previous approaches is that we use a form of the continuity 


equation that is linearized only in the velocity but fully nonperturbative in the mass density. Since reconstruction is 
most useful on small scales where perturbation theory breaks down, this nonperturbative motivation is very valuable. 

An extension that goes beyond the scope of this paper would be to use better approximations for the velocity e.g. by 
using forward-modeling techniques or iterative estimators. 


B. Halos 


Our derivation easily extends to biased tracers like halos if we assume that they define a density field that exists 
even at early times and that halos (or density peaks) are conserved over timej^j The continuity equation for the halo 
density is then 


drjSh + V • [(1 + J/tjv/j] — 0, 


(18) 


so that 


d v 8 h = - V ■ v fc - v h • - S h V ■ v h . (19) 

We now proceed to approximate the halo velocity. Assuming that the halo velocity follows the dark matter (DM) 
velocity with some non-stochastic velocity bias, = 6„v, and approximating the DM velocity as above using the 
linearized continuity equation, gives 


Vh(k) = 6„v(k) = f'Hb v l ^W R {k)5( k). 


( 20 ) 


In practice we only know the halo density Sh and not the DM density 8, but ignoring higher order bias the two can 
be related by Sh = biS so that 


h 7k 

v h (k) = fn^ ¥ w R (k)s h (k). 


( 21 ) 


The reconstruction algorithm is thus the same as for the DM case, with an additional prefactor b v /b\, which rescales 
the halo density to the halo velocity, 

b r 


’J 1 L 


Cegs( x ) = v) - A?? fU — s h (x, 77 ) ■ V8 h (x, 77 ) + 8 h (x, rj)8 h>R (x , 77 ) 


( 22 ) 


where we dropped the term that reverses linear time evolution as in Eq. ( |15[ ) 

The power spectrum of the reconstructed halo density (22) is a combination of halo 2-, 3- and 4-point statistics, 
similar to the DM case of Eq. 


( |17| ). The fiducial bias parameters used for the reconstruction enter simply as coefficients 
of halo cross-spectra. Changing these fiducial bias parameters would simply up- or down-weight the corresponding 
cross-spectra, making it computationally trivial to obtain the reconstructed power spectrum for varying fiducial 
bias parameters if all cross-spectra are pre-computed once. Such changes are more cumbersome for Lagrangian 
reconstructions, where for every new fiducial value of a bias parameter all particles need to be displaced by the 
corresponding updated displacement field. 


3 We do not discuss the validity of this assumption, but see e.g. [251 (26| for related recent discussions. We also leave the extension to 
nonlinear and nonlocal halo bias to future work, noting that we focus on dark matter in this paper. 
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III. EULERIAN F2 RECONSTRUCTION TO REMOVE 2ND ORDER PERTURBATIONS 

While the previous section was based on linearizing only the velocity but not the mass density in the continuity 
equation, we now proceed by linearizing both the mass density and velocity in all equations of motion for the DM 
fluid (i.e. the continuity, Poisson and Euler equation). Working at second order, we construct a different Eulerian 
reconstruction algorithm that aims to reverse the full second order part of the nonlinear density. As shown in 
Appendix [C] the same algorithm can be derived slightly more formally by applying the Newton-Raphson method to 
find the linear density that generates a given observed density at second order. 


A. Dark Matter 


Solving the equations of motion for DM perturbatively in the density gives the solution S = <5o + 5^ + O(Sq) where 
<5o is the linear density and the second order part is given as a sum of growth, shift and tidal term (e.g. mmmimy- 


<5 (2) ( x ) = ^o 2 (x) - ^o(x) ■ V($o(x) + ^-A'o(x) • 


growth 


shift 


tidal 


The linear displacement field in 


? k 

*o(k) = ^o(k) 


and the tidal term is given by contracting the tidal tensor 

4A) = - ij'f 


w 


to 


Ko(x) = 2 Kq ( x ) Kq 3 (x), 

where 5^ is the Kronecker delta. In Fourier space this becomes 

Kq (k) = J P 2 (k! • k 2 )<5 0 (k 1 )<5 0 (k 2 ), 
where P 2 denotes the l = 2 Legendre polynomial, 

P2(M)=4(V-yV 


The second-order density of Eq. (23) is often written as a convolution in Fourier space, 

r 7 r* 

d (2) (k) = J —^F 2 (k 1 ,k-k 1 )5 0 (k 1 )5o(k-k 1 ) = J F 2 (k 1 ,k 2 )S 0 (k 1 )S 0 (k. 


(23) 


(24) 


(25) 


(26) 


(27) 


(28) 


(29) 


where the shorthand notation f k is defined by Eq. (|3| and the symmetrized F 2 kernel contains growth, shift and tidal 
parts analogously to Eq. 


A 2 (k-j, k 2 ) = 


17 

21 


h k2 
k 2 h 


ki • k 2 + 21 P 2 (k 1 ' k 2 )- 


(30) 


To restore linear information from the nonlinear observed density we try to esti mate the second order part (23) and 
subtract it outj^] In practice we do not know the linear density do that enters Eq. (231, but we can approximate it by 


4 The Fourier convention used here differs from e.g. m so that Vl'here = — ^there- The convention used here agrees with na. 

5 More formally, one can model the nonlinear density in terms of the linear one as S = <5o + ^o] + • • • • Then the goal is to get 

<5o from S. This can be done with the Newton-Raphson root-finding method, see Appendix [C] The second step of this corresponds to 
subtracting F 2 [5, 5] from <5, which agrees with the intuitive picture of subtracting out nonlinearities from the nonlinear density. 
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the smoothed nonlinear density because linear and nonlinear densities agree on large scales. There is some freedom 
which fields in Eq. (23) should be smoothed or linearized. Motivated by the asymmetric smoothing that followed from 
the continuity equation in the last section, we choose to smooth one of the two fields entering the growth term S 2 and 
the tidal term K 2 , and we only smooth SI/ but not VS in the shift term. This defines the Eulerian F 2 reconstruction 
‘EF2’: 


<5 ef2(x) = S(x) 


^S R (x)S(x) 

growth 


s(x) • V<5(x) 

shift 



tidal 


(31) 


where Sr is the smoothed nonlinear density, s is the negative smoothed Zeldovich displacement given by Eq. (12), 
and Aj/(x) = | is a partially smoothed tidal term where AT*- 7 is defined as in Eq. (26) but using the nonlinear 

density S instead of the linear density <5o, and is defined in terms of the smoothed nonlinear density Sr instead of 

5o. 

The Eulerian F 2 reconstruction of Eq. (31) is similar to the growth-shift reconstruction of Eq. (16); the shift terms 
s • VS are equal in both methods, while the growth term SrS has a slightly different coefficient, 17/21 instead of 1, 
and the F 2 reconstruction subtracts an additional tidal part. Both methods agree with the unreconstructed density 
on very large scales because the combination of quadratic fields that is subtracted vanishes there (this follows by 
considering the squeezed limit of the corresponding kernels in Fourier space). 


B. Halos 


It is straightforward to extend this algorithm to halos if a specific bias relation between DM and halos is assumed. 
For example, if only linear and second order bias are taken into account (neglecting tidal and velocity bias), we have 


S h {x) = h [<S„(x) + F 2 [S 0 , <S 0 ](x)] + b 2 [5 0 2 (x) - (5 0 2 (x))] + 0(<5 O 3 


(32) 


where we wrote F 2 [5q, ^o]( x ) = S^ 2 \x) for the Fourier transform of Eq. (29). At linear order <5o = Sh/b \, so that 
subtracting second-order parts from the nonlinear halo density gives (ignoring smoothing, which should be included) 


^ef 2 (x) = S h (x) - ^F 2 [S h ,S h ](x) - | [S 2 (x) - (<5 2 (x))] 


(33) 


This is an estimator for the linear part biSo of the halo density given the nonlinear halo density Sh- Note that this is 
similar to the EGS reconstruction for halos given by Eq. (22), but has slightly different coefficients and an additional 
tidal term coming from the F 2 kernel. 

The algorithm can be extended to redshift space by subtracting the second order part of the redshift space density, 
which can be obtained by modifying the real space F 2 kernel to a redshift space version Z 2 (see for a review). 
The form of the resulting terms is similar to that of real space growth, shift and tidal terms from the F 2 kernel, 
but derivative operations differ in detail because of redshift space distortions. Similarly, the EGS algorithm of the 
previous section could be extended to redshift space by inverting the redshift space continuity equation. However, 
we leave these extensions for future work and restrict ourselves to reconstructions for DM in real space in this paper, 
which is somewhat cleaner and simpler. It should be noted that in contrast to our Eulerian algorithms, the standard 
Lagrangian BAO reconstruction algorithm, introduced almost 10 years ago (8|, has been extended and applied to 
redshift space, e.g. in reference lEIl- 


iv. LAGRANGIAN RECONSTRUCTIONS 
A. Possibilities 

The standard BAO reconstruction algorithm, proposed by [Sj and denoted Lagrangian growth-shift ‘LGS’ recon¬ 
struction in this paper, works as follows. First, the particles of a clustered catalog are displaced by the negative 
Zeldovich displacement s to obtain a ‘displaced’ catalog, whose density is denoted S,i [s]. Similarly, particles of a 
random catalog are shifted by the same displacement field to give a ’shifted 1 catalog with density 6 S [s]. Then the 
density difference <5d[s] — <5 s [s] gives the reconstructed density and suppresses only the nonlinear but not the linear 
part of the density (see Appendix [A| . 












As a simple extension we also consider clustered and random catalogs shifted by the positive Zeldovich displacement, 
—s. We denote the resulting densities <5d[—s] and <5 S [—s]. At first order in 5o we have (following from Eqs. (A3)-(A5) 
in Appendix |A| and ik ■ s(k) = Wn(k)5 0 (k)) 


4 1 ) [s](k) = [l-W fl (fc)],5o(k) (34) 

5^[s](k) = -W R (k)5 0 (k) (35) 

4 1 ) [-s](k) = [l + W fl (fc)],5o(k) (36) 

6^(-s](k) = W R (k)6 0 (k). (37) 


There are several possible combinations of these four shifted/displaced catalogs that change the nonlinear part of the 
density while keeping the linear part unchanged (for arbitrary smoothing): 


6l^ s = S d [ S }-6 s [s}=S 0 + O(S 2 0 ) (38) 

\ {^d[ s ] + 3d[— s]} = + O(S 2 0 ) (39) 

^d[ s] + [—s] = So + 0(5 o) (40) 

5 s [s] + [—s] = + 0(51) (41) 

Crr = <5 - C{5 S [s] + 5 S [- 8 ]} = So + 0(51) (42) 

*d[-s]-*»[-s]=«o + 0(<yg), (43) 


where c in Eq. (421 is an arbitrary constant. Beyond first order these combinations generally differ from each other. 

For reconstructions of the linear density, the nonlinear contributions to the nonlinear density should be suppressed, 
particularly the growth and shift part. Modeling all combinations by 2LPT shows that this is only the case for two 
combinations (see Appendix A4): First, the LGS reconstruction algorithm of Eq. (38 1 , which is just the standard 
BAO reconstruction algorithm. Second, a new Lagrangian random-random ‘LRR’ reconstruction algorithm, which is 
given by Eq. ([42]) if we choose c = 1/2 (see Appendix | A 4[), 


& = ^{<Us]+<y-s]}. 


(44) 


In this case, the reconstructed density is obtained as follows. A random catalog is first shifted by the negative Zeldovich 
displacement s to get 6 S [s]. The original random catalog is also shifted by the positive Zeldovich displacement —s to 
get <5 S [—s]. Then the mass densities of these two shifted random catalogs are computed and added together. Finally, 
half of that sum is subtracted from the original density of the clustered catalog to obtain the reconstructed density. 


B. Lagrangian Growth-Shift and Random-Random Reconstructions 


We now discuss the LGS and LRR algorithms in more detail. The theoretical 2LPT expression for the standard 
LGS method is (see Eq. (A32), using superscript ‘th’ to denote theoretical densities) 


Ccs th (k) = <5 th (k)- [ 

•Ik; 


f* 

fci <* 

L 

1 + — ki • k 2 

L fc 2 J 


Wtt(fc2)<5o(ki)<5o(k2), 


or more succinctly in configuration space, 

^LGS h ( x ) = <5 th ( x ) - 6o,it(x)<5 0 (x) - s 0 (x) • V5 0 (x). 


(45) 


(46) 


Here, 5q,r is the smoothed linear density and So is the negative smoothed Zeldovich displacement computed from the 
linear density do analogously to Eq. (12). Thus, in this 2LPT model, LGS reconstruction subtracts growth and shift 


terms from the nonlinear density 5, so that in 2LPT theory the reconstructed density of the Lagrangian growth-shift 
method of Eq. (38) agrees with that of the Eulerian growth-shift method of Eq. (16), derived from the continuity 


equation for the nonlinear mass density. 

Note that the statement here is that the reconstructed densities agree at second order LPT. For the reconstructed 
power spectrum this means that the tree level and ‘ 22 ’ loop contributions agree between the methods, whereas the 
‘13’ loop correction (which is of the same order as the ‘22’ contribution) can differ between these algorithms because 
the reconstructed densities may differ at third order, which we do not model here; see Section |VII A 3| for further 
discussion of this point using simulations. 
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Details of the 2LPT calculation that leads to Eqs. (451 and (46) are given in Appendix [A| These results only 
hold after introducing a new correction term that arises from modeling the displaced clustered catalog such that the 
displacement field s is evaluated at Eulerian instead of Lagrangian positions; see Appendices [A] and [B| 

Similarly, modeling the random-random reconstruction of Eq. (44) with 2LPT gives (see Appendix |A 4|) 


rrec,th 

°LRR 


(k) = ^(k)- 


'k; L 


2 1 (k x k 2 \s ~ l n ~ ' 

2 + 2 ( + JT^j ki • k 2 + -P 2 (ki • k 2 ) 


Wfl(fci)tE fl (fc 2 )J 0 (ki)(5o(k 2 ), (47) 


The square brackets correspond to the part of the F 2 kernel that is generated by [l/ 1 )] 2 in 2LPT; see Eq. (A23). In 
configuration space this becomes 


CRR h W = <5 th M - ^o.r(x) - s 0 (x) . V5 o,r(x) - ^A' 2 m (x), 


(48) 


where Kg RR (x) = § K l g R Kg R is defined in terms of the smoothed linear density. 


' 2 fx'l - 2 
L 0 ,RB\ 1 ~~ 2 

The nonlinear density before reconstruction can be modeled by 

^ ■ 17 


<5 (x) = (S 0 (x) + — <5 0 (x) - (x) • Vd 0 (x) + —K 0 (pc) 


(49) 


at second order LPT or SPT. On large scales, <5o ~ ^o,r ~ and \&o ss —s 0 « — s. Thus, both reconstruction 
methods LGS and LRR cancel the shift term • V<5 exactly on large scales as desired. 

However, even at second order the reconstructions do not perfectly recover the linear density <5o - For the LGS 
method, the residual second order part of the density after reconstruction is 


Cgs (k) - <S 0 (k) = 


17 

21 


- W R (k 2 ) 


+ [1 — W R {k 2 )\ Tykq • k 2 + — P 2 (ki • k 2 ) ldo(ki)<5o(k 2 )- (50) 

k 2 21 


On large scales, where W R {k\) — > 1, reconstruction reduces the growth part from 17/21 to —4/21 and it fully removes 
the shift term, but it does not change the tidal part of the nonlinear density. On smaller scales, where HOi/fci) —> 0, 
reconstruction does not change the second order part of the density, which is reasonable because the linear displacement 
field cannot be reliably estimated from the nonlinear density any more on these scales, and still attempting to do so 
may inadvertently enhance nonlinearities j^] 

For the LRR method the residual second order density is similar, 


O h (k)-<5 0 (k) = 


- ^w R (h)w R (k 2 ) 


+ [1 — W R {k\)W R (k 2 )} — 1 -ki • k 2 
& 2 


^ - ^W fl (fci)W fl (fc 2 ) 


P 2 (kr-k 2 ) K(ki)5 0 (k 2 ). 


(51) 


For large scale modes, the growth term is reduced from 17/21 to 3/21, the shift term is removed entirely, and the tidal 
term is changed from 4/21 to —3/21; on small scales the second order part is again unchanged. Note that the LRR 
method involves two smoothing kernel factors that smoothen both displacement and density fields. This aggressive 
smoothing may remove small-scale density modes that could actually still be used if they were displaced by large-scale 
displacements, at least for applications like measuring the BAO scale. 

On the very largest scales, k —> 0, the second order residual densities given by the right hand sides of Eqs. (50) and 
(51) vanish. To see this, note that for k <C ki,k 2 , we have k 2 = k ki « —ki corresponding to a squeezed triangle 
configuration, so that k <C fci ~ k 2 . Then we get P 2 (—1) = 1. For ki,k 2 < i? -1 , W R (ki) = 1, and the right hand 
sides of Eqs. ([50|) and (51) become —4/21 + 4/21 = 0 and 3/21 — 3/21 = 0, respectively. In the regime ki,k 2 > Rr 1 , 


W R (ki) = 0, and the right hand sides of Eqs. (50) and (51) both become 17/21 — 1 + 4/21 = 0. 


Reversely, this fact could be exploited by applying reconstruction only to modes or local environments for which the velocity or 
displacement can be estimated reliably eh. 
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C. Lagrangian F2 Reconstruction 


The 2LPT theory expressions (|45|) and (47) for the LGS and LRR methods (or the residual equations (50) and 


(51)) can be used to form a linear combination of these two reconstruction algorithms that removes the full F 2 kernel 
(oil scales that are not smoothed out) and still leaves the first order part unchanged^ 


<Tf 2 = ^LGS + PlKK = ^ + ^[s] - ^S t [s) 
where ‘LF2’ stands for Lagrangian F 2 reconstruction. Indeed, at second order in LPT 


(52) 


C F c f(k) = nk)^ / Wn(k 2 ) 




3 4 

7 + jW R (k i) 


+ 7^W R (ki)P 2 (ki ■ k 2 ) ^o(ki)<5o(k 2 ), 


l“ kl 

k 2 


(53) 


where the kernel in curly brackets equals F 2 for large scales ki < R~ x , where W R (kf) = 1, i.e. the full second order 
F 2 part of the nonlinear density is subtracted on large scales as desired. In 2LPT this Lagrangian F 2 algorithm is 
equivalent to the Eulerian F 2 method (at the level of reconstructed densities, and up to the freedom which fields to 
smooth). 


V. OVERVIEW OF RECONSTRUCTION ALGORITHMS 

Before discussing simulations and numerical results, we briefly summarize all reconstruction algorithms in Table [I] 
and in the following list. The Eulerian algorithms that operate directly at the field level are: 

1. EGS: Eulerian Growth-Shift Reconstruction (Section\U) 

Nonlinear growth and shift parts are subtracted from t 


le unreconstructed density <5, 
<*egs( x ) = 5 ( x ) - M X )<K X ) ^ s ( x ) ' V6(x), 


(54) 


where S R is the smoothed density and s is defined in Eq. (12). This follows from linearizing the velocity but 


not the mass density in the nonlinear continuity equation. The reconstructed power spectrum combines the 
unreconstructed power and specific 3- and 4-point statistics of the unreconstructed density, see Eq. ©• 


2. EF2: Eulerian F2 Reconstruction (Section III) 

Aiming to reverse the full second-order S ^ 2 ' ~ F 2 5 * S part of the unreconstructed density, we form 

^ef 2 ( x ) = <5( x ) ^ ^M X M( X ) - s(x) • V(5(x) - ^-/v^(x). 


(55) 


This is similar to the growth-shift reconstruction (|54|) but an additional tidal part is subtracted and the 


growth term has a slightly different coefficient. The reconstructed power spectrum is again a combination of 
specific 2-, 3- and 4-point statistics of the unreconstructed density similar to Eq. ©• 

3. ERR: Eulerian random-random, reconstruction 

This Eulerian algorithm resembles the Lagrangian random-random LRR algorithm (see Eq. ©), 


^ERR = <K X ) - 3<5 r( x ) - s(x) • VS R (x) - -I< RR (x) 


(56) 


Additionally, we consider the corresponding Lagrangian reconstruction algorithms that require displacing objects 
in catalogs: 


7 Ignoring smoothing, there is exactly one linear combination ci<5£q S + C 2 ^£r R that works. To keep the correct shift term we need 
ci -f- C 2 = 1. The growth term imposes ci + |c 2 = |y. The tidal term is fixed by the fact that the kernel needs to vanish in the squeezed 
limit (or it imposes ^C 2 = ^y). The only solution is ci = y and C 2 = j. 
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Eulerian 

Lagrangian 

Growth-Shift 

£egs( x ) = <H X ) - Sr(x)5(k) - s(x) ■ VS(x) 

Cos = Ms] - <Ms] 

Full F2 

^EF 2 ( x ) = <K X ) - - s(x)V<5(x) - A -K%(x) 

rrec 3 rrec 1 4 rrec 

°LF2 — 7°LGS ~v y°LRR 

Random-Random 

& = <5( x ) - §4( x ) - s(x) ■ VMx) - \K 2 rr (x) 

Crr = <5-! (<Ms]+<M-s]} 


TABLE I. Overview of the reconstruction algorithms considered in this paper. Lagrangian reconstructions involve displacements 
of objects which is not the case for Eulerian reconstructions. Methods in the same row agree at second order LPT. The standard 
BAO reconstruction algorithm is called Lagrangian growth-shift reconstruction and listed on the upper right. 


1. LGS: Lagrangian Growth-Shift Reconstruction (Section \IV B ) 

This is the standard BAO reconstruction algorithm proposed in [Sj: Objects in clustered and random catalogs are 
displaced by the large-scale negative Zeldovich displacement field s, and the difference between the corresponding 
‘displaced’ and ‘shifted ’densities (5d[s] and <5 s [s] gives the reconstructed density, 

Cgs = Ms] - W (57) 


2. LF2: Lagrangian F2 Reconstruction (Section IVC) 

The LGS and LRR methods are combined such that the second order F 2 part of the unreconstructed density is 
subtracted on large scales, 


3 4 

rrec _ rrec i rrec 

°LF2 = j°LGS “T y°LRR- 


(58) 


3. LRR: Lagrangian Random-Random reconstruction (Section IVB) 

This involves two oppositely shifted randoms: Randoms are shifted by the negative Zeldovich displacement, <S s [s], 
and by the positive Zeldovich displacement, 5 S [—s], and their mean is subtracted from the observed density, 


S% R = 8--{5 a [B] + 8s[-B]}- 


(59) 


The result agrees with the linear density at first order and suppresses the nonlinear growth and shift term at 
second order. 


As explained before, the 2LPT model for the reconstructed density of each Eulerian algorithm agrees with the 
corresponding Lagrangian-reconstructed density, although the algorithms are operationally very different. All recon¬ 
struction algorithms leave the power spectrum unchanged on the very largest scales and reverse nonlinearities on 
smaller scales in different ways. 


VI. SIMULATIONS 

To test the performance of the above reconstruction algorithms, we ran several N-body simulations with the 
FastPM code [30], derived from the parallel COLA m implementation of [nng. This is a particle-mesh (PM) 
code that uses linear time stepping in the scaling factor a to produce accurate large scale structure at a fraction of 
the total computing time of a typical TreePM N-body simulation. We set up 2LPT initial conditions for 2048 3 DM 
particles in a box of length L = 1380Mpc//i per side and evolve them from the initial redshift z t = 99 to present 
time with 80 time steps. From z = 99 to z = 1 we use a 4096 3 PM grid; at z < 1 we switch to a 6144 3 PM grid. 
We only use a single snapshot at z = 0.55, which is taken after 52 time steps. We run three independent realizations 
by generating Gaussian random fields following an initial DM power spectrum. For each realization we also run a 
nowiggle simulation that has the same phases but is obtained from a nowiggle power spectrum where BAO wiggles are 
smoothed out. The cosmic variance caused by broadband fluctuations can thus be cancelled when comparing wiggle 
and nowiggle simulations [33] [34]. We assume a flat ACDM cosmology with fl m = 0.272, h = 0.702, erg = 0.807 for 
the N-body simulations. The convergence of the simulations is discussed in Appendix |D 1| 

For all plots that do not require wiggle and nowiggle simulations, we instead use ‘RunPB’ N-body simulations 
produced by Martin White with the TreePM N-body code of [35]. These simulations were also used in [Ml 1M1 [33- 
They have 10 realizations of 2048 3 DM particles in a box with side length L = 1380Mpc//i, and were started with 2LPT 
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initial conditions at z t = 75 for an initial power spectrum with BAO wiggles (these simulations were not run for an 
initial nowiggle power spectrum). The cosmology is flat ACDM with f \h 2 = 0.022, fl m h 2 = 0.139, n s = 0.965, h = 0.69 
and erg = 0.82, and we again only use one snapshot at z = 0.55. 

To speed up reconstructions we work with 1% DM subsamples0of clustered catalogs, where the same random subset 
of particles is selected from wiggle and nowiggle simulations. Random catalogs are generated by distributing equally 
many particles randomly in the box. Densities are computed on 512 3 grids using the cloud-in-cell (CIC) scheme. 
The CIC kernel is deconvolved from the density grids. For densities that need to be smoothed when entering in a 
reconstruction method we employ Gaussian smoothing with kernel Wn(k) = exp(— ^k 2 R 2 ). Following [TUI ITT] we 
choose R = 15Mpc//i as the fiducial smoothing scale (see Appendix D2 for different smoothing scales). 

For Eulerian reconstructions, we calculate quadratic fields like \&(x) • V<5(x) by Fourier-transforming <5(x) to k 
space, multiply the result by appropriate filters like k /k 2 , Fourier-transform the result back to x space and multiply 
fields there. These fields are then combined according to the particular reconstruction algorithm under consideration. 
Finally, we measure the power spectrum of the reconstructed density. For more details we refer the reader to [ 21] 
where spectra of the same quadratic fields were measured. 

For Lagrangian reconstructions, the displacement field s is computed in k space using Eq. (12) (see also 01]). This 
is then Fourier-transformed to x space and linearly interpolated to the locations of particles in clustered and random 
catalogs. The particles are displaced according to the interpolated displacement field and the density of the displaced 
catalogs is obtained with the CIC method. Finally we form the power spectrum of the reconstructed density. 

To cancel the broadband cosmic variance between wiggle and nowiggle simulations when post-processing measured 
spectra for plots, we compute the required combination of spectra (e.g. the fractional difference between wiggle and 
nowiggle simulations) for each realization individually, and only average the final result over realizations. Error bars 
are standard errors of the mean obtained in this way, and estimated from the scatter between realizations. Due to 
the cosmic variance cancellation these errors are typically very small and we will therefore omit them in most plots 
for clarity. 


VII. NUMERICAL RESULTS 

We now discuss how well the reconstruction algorithms perform in simulations. For clarity we first discuss the 
performance of the Eulerian and Lagrangian growth-shift reconstructions EGS and LGS in isolation and study in 
detail where additional BAO information comes from. This is then followed by a more comprehensive comparison 
of all algorithms in Section |VII B| To be quantitative, we only consider power spectra here. Appendix [B] compares 
the algorithms at a somewhat more qualitative level using 2D density slices, showing that the density changes due to 
Eulerian and Lagrangian reconstructions are very similar and agree with 2LPT expectations at the field level. 


A. Growth-Shift Reconstructions 

1. Overall Comparison 

Fig-0 shows power spectra before and after reconstruction. For the Lagrangian LGS method, the broadband 
shape is similar before and after reconstruction, whereas the Eulerian EGS reconstruction adds additional broadband 
power at high k. This is only true for the particular smoothing that we use for the EGS method. E.g. if we used 
symmetric smoothing where all fields entering into quadratic fields are smoothed rather than just one of them, then 
the broadband of the reconstructed power agrees with that before reconstruction. However, symmetric smoothing 
yields significantly lower BAO signal-to-noise than asymmetric smoothing, which makes sense because asymmetric 
smoothing allows small-scale fluctuations to be shifted by large-scale flows, whereas symmetric smoothing erases all 
small-scale information. 

To highlight the BAO signal per fc-mode, Fig. [2] shows the fractional difference of spectra between simulations 
with and without wiggles. After reconstruction, this fractional difference has more pronounced wiggles than before 
reconstruction and is closer to linear theory, implying an enhanced BAO signal. The enhancement is similar for 
Eulerian and Lagrangian reconstructions. Most (but not all) of the signal can also be obtained by forming the 
cross-spectrum between reconstructed and unreconstructed densities, as shown by non-filled symbols in Fig. [2] 


For R = 15Mpc /h smoothing, the relative performance of the reconstructions does not change qualitatively for 0.1% and 0.01% clustered 
subsamples, if the random catalogs for Lagrangian reconstructions contain as many particles as a 1% clustered subsample. 
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k [/i/Mpc] 


FIG. 1. Total DM power spectra before reconstruction (black squares), after standard Lagrangian growth-shift LGS recon¬ 
struction by particle displacements (circles), and after Eulerian growth-shift EGS reconstruction based on combining 2-, 3- and 
4-point information (stars), for 1% subsamples of simulations with BAO wiggles. Filled symbols show density auto-spectra, 
while open symbols show cross-spectra between the densities before and after reconstruction. Results are averaged over 3 
realizations and error bars correspond to the standard error of the mean (they are smaller than the symbols except at very low 
k where their size is similar to that of the symbols). 



k [h/ Mpc] 


FIG. 2. Illustration of the BAO signal: Fractional difference of the power spectra in Fig. [I] between wiggle- and no-wiggle 
simulations. Again, error bars are smaller than the symbols. 


Fig-! shows the signal-to-noise-squared per k mode estimated by assuming that the BAO signal is given by the 
difference of power spectra between wiggle and nowiggle simulations, and approximating the covariance with the 
Gaussian expectation 

cov(P(k), P(k')) = S ktk i— - j—r(P(k)) 2 , (60) 

where the number of modes per fc-bin and the power spectrum expectation value are computed from the simulations; 
see Appendix |D 3| for justification of this approximation. Fig. [3] shows that the BAO signal-to-noise-squared is clearly 
enhanced after EGS and LGS reconstructions. The algorithms perform equally well for k < 0.15ft./Mpc. On smaller 
scales, the Lagrangian LGS algorithm performs slightly better than the Eulerian EGS algorithm. This is caused by 
the larger broadband noise of the EGS algorithm in this regime (see Fig. [I]). 
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CNI 




FIG. 3. Signal-to-noise-squared for the BAO wiggles as a function of wavenumber k. The signal is given by the difference 
of the spectra specified in the legend between wiggle and nowiggle simulations, S = Pwiggie — Pnowiggie- The squared noise is 
approximated by Eq. (601. Error bars are not shown for clarity. 



To better quantify the overall performance, Fig. [4] shows the cumulative BAO signal-to-noise-squared as a function 
of fc max . The Lagrangian LGS reconstruction improves the BAO signal-to-noise by a factor 1.4 for fc max ~ 0.4/i/Mpc. 
The Eulerian EGS algorithm yields \/l84/202 ~ 95% of the BAO signal-to-noise of the LGS method. Thus, Eulerian 
and Lagrangian growth-shift reconstructions perform very similarly, differing by less than 5% in their total BAO 
signal-to-noise (see also Table [TT] later). 

Using cross-spectra between reconstructed and unreconstructed densities instead of auto-spectra of reconstructed 
densities degrades the cumulative signal-to-noise by ~ 5% for each of the two reconstruction methods. 


2. 3-point vs point 


Does most of the additional BAO information gained by reconstruction come from the 3-point or the 4-point 
function of the unreconstructed density? To answer this, we write the reconstructed density as 

<5 rec = J + ( 6 iec - 6), 


(61) 
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k [/i/Mpc] 


FIG. 5. Illustration of 3- and 4-point contributions to reconstructed power spectra. The full reconstructed power spectra (circles 
and stars) can be obtained by adding the 3-point part (thin solid) and 4-point part (dash-dotted) to the pre-reconstruction 
power spectrum (black squares); see Section VIIA 2 for details. The curves show differences of the spectra specified in the 
legend between wiggle and nowiggle simulations, divided by the theoretical linear no-wiggle density power spectrum. At all k 
error bars estimated from the scatter between the 3 realizations are roughly the same size as the thickness of the colored curves 
and smaller than the size of the circle, star and square symbols; they are not shown for clarity. 


where <5 is the unreconstructed density. For Eulerian reconstructions, the density change due to reconstruction is 
quadratic in the unreconstructed density, ( S rec — S) ~ 0(6 2 ). The reconstructed power spectrum therefore splits into 
2-, 3- and 4-point parts as 

(<5 rec |<5 rec ) = (5\5) + 2((<T ec - <5)|J) + ((4 rec - S) \ (6 iec -6)) . (62) 

2pt 3pt 4pt 


For example, the cross-spectrunr between (i5 rec — S) and S is a cross-spectrum between a field quadratic in S and S 
itself, and is therefore given by a particularly integrated bispectrum of <5 [24] . 

Simulation measurements of the terms on the right hand side of Eq. (62) are shown in Fig. [5 For EGS reconstruction, 
the enhanced BAO signal is almost entirely due to the 3-point part of the signal (thin solid), whereas the additional 
signal from the 4-point (dash-dotted) is almost negligibly small. 

For Lagrangian reconstructions the difference between reconstructed and unreconstructed densities is in general 
not simply quadratic in the observed density any more, because the reconstructed density is obtained by calculating 
new densities from shifted and displaced catalogs that cannot easily be related to the density of the original catalog. 
It is therefore less transparent to study the role of 3- and 4-point signal in this case. Still, modeling Lagrangian 
reconstruction from a Lagrangian theory perspective shows that reconstructed and unreconstructed densities agree 
at first order in the linear density, so that the leading-order difference between reconstructed and unreconstructed 
densities is quadratic in the linear density; see Eq. (46) and Appendix[Aj The spectra on the right hand side of Eq. (62) 
are therefore still related to 3- and 4-point contributions if third order corrections are ignored. Fig. [5] shows that for 
LGS reconstruction these spectra are similar to EGS reconstruction, indicating that the enhanced BAO signal from 
LGS reconstruction is also mostly due to the effectively included 3-point signal. 


3. 13 and 22 contributions to the 3-point part of the reconstructed power spectrum 


To further elucidate the origin of the enhanced BAO signal after reconstruction we split the observed density 5 into 
linear and nonlinear part 


4 = S 0 + (5 - < 5 0 ), 


( 63 ) 
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k [/i/Mpc] 


FIG. 6. Separation of the reconstruction 3-point part 2((5 rec — <5)|<5) (thin solid) from Fig. [5] into ‘13’ contribution (dashed) 
and ‘22’ contribution (dotted) according to Eq. (641, measured from simulations for standard LGS reconstruction (yellow) and 
EGS reconstruction (green). The linear density <5o in the cross-spectra is approximated by the density of the simulation initial 
conditions, rescaled to z = 0.55 with the linear growth factor. Adding dotted to dashed curves gives the thin solid curve of 
the same color. The curves show differences of the spectra specified in the legend between wiggle and nowiggle simulations, 
divided by the theoretical linear no-wiggle density power spectrum. Error bars of all colored curves estimated from the scatter 
between the 3 realizations are roughly the same size as the thickness of the curves at all k ; they are not shown for clarity. 



0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 

k [/i/Mpc] 


FIG. 7. Same plot as before but including broadband information by plotting spectra specified in the legend in wiggle simulations 
divided by the theoretical linear no-wiggle density power spectrum. Error bars are roughly the same size as the thickness of 
the curves at k < 0.15/i/Mpc and twice that size at higher k , but they are not shown for clarity. The plot shows the average 
over the 10 RunPB N-body realizations available to us. For the 3 FastPM realizations that we ran for the previous plots this 
looks very similar but slightly noisier because of the smaller number of realizations. 
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where <5o is the linear density. Then, the 3-point-like term from the last section splits into 

((<5 rec - $)!<*> = ((<5 rec - <5)|<5 0 ) + ((S rec -6)\(5- S 0 )). 


^A(3)5 0 )+e>(5g) 


-(Al 2 )5( 2 ))+0(5g) 


Expanding 


i5 rec - 5 = A (2) + A (3 > + ■ 


(64) 


(65) 


where A^") is n-th order in the linear Gaussian density fig, implies that the first term on the right hand side of Eq. (64) 
corresponds to a contraction of with <5o at leading order in So', we call this the ‘13’ contribution. The second 
term on the right hand side of Eq. (641 corresponds to a ‘22’ contribution at leading order. 


To see this more explicitly, we expand the field change due to reconstruction up to third order in <5o, 

<5 rec - <5 = A (2) + A (3) ( 66 ) 

= I (2^) 3 fo(k-k 1 -k 2 )^ s) (k 1 ,k 2 ) ( 5o(ki) l 5o(k 2 ) 


'ki,k 2 


/k | ,k 2 .k 3 

imetric 1 

consideration. Then, the 13 contribution to Eq. (|64| becomes 


(2<) 3 Mk - kr - k 2 - k 3 )Dcj s) (k 1 , k 2 ,k3)<5o(ki)^(k 2 )^(k 3 ), 


(67) 

( 68 ) 


(s) (s) 

where H 2 and D 3 ' are symmetric kernels whose explicit form depends on the reconstruction algorithm under 


(69) 


(A( 3 )(k)5 0 (k')> = 3(27r) 3 Mk + k')Pii„(£0 [ D^ ) (k 1 ,-k 1 ,k)P lin (fc 1 ), 

J ki 

and similarly the 22 contribution is 


(A< 2 >(k)<f< 2 >(k')) = 2 (27r) 3 <5o(k + k') 


(27r) 3 fe(k-k! -k 2 )P^ ) (k 1 ,k 2 )P 2 (k 1 ,k 2 )P lin (fc 1 )P lin (fc 2 ). (70) 


'ki,k 2 


Figures [ 6 ] and [7] show the spectra on the right hand side of Eq. (64) measured from simulations. The linear density 
<5o is approximated by estimating the density of the initial conditions of the simulations and rescaling it to z = 0.55 
with the linear growth factor. The measu red 13 term has pronounced BAO wiggles, which can be attributed to the 
fact that the theoretical 13 term of Eq. (69) is proportional to Pu n (fc) which has BAO wiggles]^] In contrast, the 
measured 22 contribution in Fig. [ 6 ] oscillates out of phase with the linear BAO wiggles, which shifts the BAO scale 
inferred from scales k < 0.15/i/Mpc. This can be understood perturbatively from the form of the theoretical 22 
expression in Eq. © imnziEsi. At higher k, the 22 contribution is opposite in phase to the linear BAO wiggles 
and thus reduces the enhancement of BAO wiggles from the 13 contribution j 33 ] 

In summary, most of the additional BAO signal from reconstruction is due to the 13-type part of the 3-point 
contribution to the reconstructed power spectrum, whereas the 22-type 3-point and the 4-point do not enhance the 
signal much. 

Figures [ 6 ] and [ 7 ] also compare the 13 and 22 contributions between the Eulerian and Lagrangian growth-shift 
reconstruction algorithms. These algorithms have the same 22 contributions at k < 0.15/i/Mpc, which is expected 
because they are equivalent when modeled with 2LPT (see Appendix 0 . In contrast, the 13 contributions differ 
even at low k, which indicates that the algorithms differ at third order in the linear density. This is not surprising 
since they involve rather different transformations of the unreconstructed density. Since the 13 contribution is most 
important for improving the BAO signal it may be worth considering third order reconstruction algorithms that 
subtract F :i 5 * S * S correctly. However, this goes beyond the scope of this paper, and it should be noted that for 
Eulerian reconstructions such an extension might suffer from shot noise issues because e.g. 5 3 (x) is very sensitive to 
the largest peaks of the density. In this case, the standard reconstruction method may be a more elegant solution, 
because it successfully reverses nonlinear BAO smoothing beyond second order without the need to cube any Eulerian 
fields. 

All spectra in Fig. [ 7 ] vanish on very large scales (k < 0.05/i/Mpc) as expected because both reconstructed and 
unreconstructed densities approach the linear density in this regime. 


9 This makes the approximation that the convolution of P\in with Ds in Eq. ( |69| ) yields a broadband power that affects the wig gles only in 
a sub-dominant way. Also, note that even at sixth or higher order in 5o, the do(k') factor on the l.h.s. implies that Eq. \69\ will always 
be of the form Pu n (k) times some integral over weighted power spectra, which can lead to further enhancements of the BAO wiggles. 

10 If the kernels were known, BAO information from the convoluted 22 contributions could still be extracted, but in practice the kernels 
are only known in perturbation theory which starts to break down on scales where reconstruction helps most, so that more information 
can be extracted by the standard approach of fitting an oscillatory function on top of a broadband shape out to smaller scales. 
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FIG. 8. Separation of EGS reconstruction in contributions from growth term S 2 (blue) and shift term s- V<5 (red), see Eq. (171. 
The plot shows the cumulative BAO signal-to-noise-squared of spectra specified by the legend. Green stars show full EGS 
reconstruction, colored solid show reconstructions if only the growth (blue) or only the shift term (red) is used. Dashed lines 
show the signal-to-noise-squared of individual 3-point contributions to the full EGS reconstruction. Dash-dotted and dotted 
show 4-point contributions. 


4- Growth vs Shift Term 


For Eulerian growth-shift reconstruction, = S — S 2 — s- VS, the power spectrum of the reconstructed density can 

be split into six cross-spectra between unreconstructed density <5, growth term S 2 and shift term s • V<5, see Eq. (171. 
Fig. [8] shows the cumulative BAO signal-to-noise from each of these contributions^] The cross-spectrum between 
shift term and the density is most important, which makes sense intuitively because this term corresponds to shifting 
back large-scale flows. If the reconstruction only included the shift term, <5 rec = <5 — s • VS (red solid in Fig. [8]), then 
the improvement over performing no reconstruction (black squares) is however only half of the improvement obtained 
from the full EGS reconstruction (green stars). Thus, the growth term is subdominant but adds very significant 
signal-to-noise so that it should not be omitted. Also, if the growth term was not included, the reconstructed power 
spectrum would not yield the linear power spectrum on large scales, so that such a shift-only reconstruction scheme 
should not be taken seriously in any case. 


B. All Reconstructions 

So far we only discussed numerical results for the growth-shift reconstructions. We will now turn to the F2 and 
random-random reconstructions introduced in Sections |III| and |IV| and compare all reconstruction algorithms. 

All six reconstruction methods are compared in Fig. 19] in terms of their cumulative BAO signal-to-noise-squared 
as a function of fc max . Table [TT] shows the total BAO signal-to-noise up to fc max = 0.4/i/Mpc. Overall, the standard 
Lagrangian growth-shift (LGS) algorithm works best. Compared against that, the Eulerian growth-shift (EGS) and 
Eulerian F2 (EF2) methods perform almost as well, missing only 5% and 6% of the total BAO signal-to-noise, 
respectively. The Lagrangian F2 (LF2) method turns out slightly worse, missing 11% of the LGS signal-to-noise. 
Both random-random reconstructions perform significantly worse, missing more than 20% of the LGS signal-to-noise. 
Note that the growth-shift and F2 algorithms perform equally well on large and intermediate scales, k < 0.15/i/Mpc, 
as expected from 2LPT modeling, and they only begin to differ on smaller scales. Thus, EGS, EF2 and LF2 only 
perform slightly worse than LGS because they partially miss highly nonlinear BAO information that the LGS algorithm 
restores slightly better. 


11 The individual signal-to-noise-squared curves in Fig.[8]cannot easily be combined by eye because there are nontrivial correlations between 
the spectra. These are taken into account for the full reconstructed signal-to-noise-squared because spectra are combined realization 
by realization (which is equivalent to combining S and the growth and shift term at the field level for each realization and then taking 
spectra of the resulting field). 
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FIG. 9. Cumulative BAO signal-to-noise-squared for reconstructed power spectra as a function of fc max measured from sim¬ 
ulations for three Eulerian reconstruction algorithms (dashed colored) and the corresponding three Eulerian reconstruction 
algorithms (solid colored). See Section [V] for a summary of the algorithms. The signal-to-noise-squared of the density before 
reconstruction (black) and of the linear density (gray) are included for comparison. 


The EGS and EF2 methods may be useful for applications because they perform almost as well as the standard LGS 
algorithm and have potential advantages due to their Eulerian nature (e.g. they can be expressed as combinations of 
2-, 3- and 4-point functions of the unreconstructed density). The LF2 method and the random-random methods are 
not competitive in terms of signal-to-noise so they should not be used. The performance of ERR can be improved 
by using asymmetric smoothing where only one of the fields entering quadratic fields is smoothed; in this case ERR 
performs almost as well as EF2 but still slightly worse. For the LRR method it is less obvious how to reduce the 
smoothing to improve it, and we do not investigate this further. 

So far we focused on the BAO signal-to-noise to assess reconstruction performance. Another important aspect of 
the standard LGS reconstruction is that it removes a shift of the BAO scale induced by nonlinear clustering 

DSHD]. 

We investigate this effect for our new algorithms as follows. For every estimated power spectrum we construct the 
fractional difference (P w iggie — P n owiggie)/Pnowiggie shown in Fig. [2] Then we fit this with the damped shifted linear 
prediction 0(k/a) exp[-fc 2 E 2 /2] for varying a and E, where 0(k) = [Peggie( fc )- P nowiggie( fc )]/ P no n wig g ie( fc ) corresponds 
to the gray line in Fig. [2j Focusing on the BAO wiggles, we fit over the range 0.083 h/Mpc < k < 0.3/z/Mpc. Since 
the sample standard deviations estimated from the scatter between three realizations are somewhat uncertain for 
individual k bins, and since errors have no broad-band scale dependence because broad-band cosmic variance is 
cancelled, we replace sample errors by their average over all k bins, which gives a constant error of 2 x 10” 4 . 

Marginalizing over the damping parameter E, the best-fit values for the BAO shift a — 1 after LGS, EGS, LF2, 
EF2, LRR and ERR reconstructions are 0.003%, 0.01%, 0.06%, 0.03%, 0.13% and 0.13%, respectively, and 0.25% 
without reconstruction. The la fitting uncertainty is roughly 0.01%. This shows that mitigation of the BAO shift 
is more efficient for reconstructions that also yield more BAO signal-to-noise. In particular, for LGS, EGS and LF2 
reconstructions the residual BAO shift is at most 0.03%, which is negligibly small for upcoming experiments. Within 
the fitting uncertainty, the LGS and EGS reconstruction do not show evidence for any residual BAO shift. 

Our results are based on simulations and it is not clear how to explain them within the framework of second order 
perturbation theory where all methods agree up to smoothing and should thus perform similarly well. The difference 
in performance could be attributed to different smoothing operations, e.g. theory expressions for the random-random 
methods involve only the smoothed density, whereas theory expressions for the other reconstruction methods involve 
also unsmoothed fields so that e.g. small-scale modes are shifted by large-scale smoothed displacements, which can 
yield improved reconstructions. However this does not explain the (small) differences between the LGS, EGS, LF2 and 
EF2 algorithms. These could be due to higher order corrections in perturbation theory, which we have not included 
in our modeling where the density was truncated at second order in the linear density. The fact that the 13 part of 
the 3-point contribution to the reconstructed power spectrum seems to be most important in simulations (see Fig. [6]) 
suggests that third-order corrections could be responsible for the slightly different performances of the LGS, EGS, 
LF2 and EF2 algorithms. We do not investigate the performance differences of these algorithm in more detail here 
because the LGS, EGS and EF2 algorithms perform almost equally well. 
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Growth-Shift F2 reconstruction Random-Random 


Reconstruction method 

LGS 

EGS 

LF2 

EF2 

LRR 

ERR 

Perfect 

NoRec 

BAO signal-to-noise 

14.2 

13.6 

12.7 

13.3 

11.1 

11.3 

17.0 

10.3 

Compared against LGS 

±0% 

-4.7% 

-11% 

-6.3 % 

-22% 

-21% 

+19% 

-27% 

Compared against NoRec 

+38% 

+31% 

+23% 

+29% 

+6.9 % 

+8.8% 

+64% 

±0% 


TABLE II. Total BAO signal-to-noise for fc max = 0.4/i/Mpc for various reconstruction algorithms (obtained from Fig. [ 9 ] based 
on simulations). ‘Perfect’ refers to the BAO signal-to-noise of the linear density and ‘NoRec’ to the BAO signal-to-noise of 
the measured nonlinear density without performing any reconstruction. The second-to-last row shows how much of the signal- 
to-noise is lost compared to performing the standard LGS reconstruction. The bottom row shows how much signal-to-noise is 
gained by reconstructions compared to performing no reconstruction. 


VIII. CONCLUSIONS 

Density reconstruction plays an important role in modern cosmology by significantly improving BAO information 
from galaxy surveys. In this paper we present three new Eulerian reconstruction algorithms that operate directly on 
the pre-reconstruction density field, and two new Lagrangian reconstructions that displace objects in catalogs. These 
algorithms are physically and operationally different from the conventionally used Lagrangian growth-shift (LGS) 
reconstruction. 

The Eulerian growth-shift (EGS) reconstruction algorithm is derived by reversing the nonlinear continuity equation, 
keeping the mass density fully nonperturbative and only approximating the velocity density in terms of its linear 
relation to the mass density. The resulting reconstructed power spectrum is given by a particular combination of 
simple 2-, 3- and 4-point statistics of the unreconstructed density. In our simulation setup, this combination yields 
the same BAO signal-to-noise as the standard LGS algorithm up to fc max = 0.15/i/Mpc, and 95% of the signal-to-noise 
of the standard LGS algorithm up to £; max = 0.4/i/Mpc. The EGS and LGS reconstruction algorithms are thus very 
similar in their ability to extract additional BAO information. 

With the goal of removing the entire second-order part of the nonlinear density, i.e. nonlinear growth, shift and 
tidal terms, we introduced Eulerian and Lagrangian F2 algorithms. The Eulerian version, EF2, performs almost as 
well as the growth-shift algorithms, missing only 6% of the BAO signal-to-noise of the standard algorithm in our 
simulations. The corresponding Lagrangian algorithm, LF2, performs somewhat worse, likely because small-scale 
modes are smoothed out too aggressively in this case. Motivated by 2LPT results, we also considered random-random 
reconstructions, where random catalogs are shifted by the positive and negative Zeldovich displacement. They turn 
out to perform significantly worse than the other algorithms. In summary, the standard LGS algorithm performs 
best, but two new Eulerian reconstruction algorithms, EGS and EF2, perform almost equally well; see Fig. [9j 

For all our Eulerian algorithms, the reconstructed power spectrum is obtained by adding to the unreconstructed 
power spectrum simple model-independent 3-point and 4-point statistics of the unreconstructed density, like (<5 2 |<5) 
or (<5 2 1 <5 2 ). This makes it very transparent how these Eulerian reconstructions yield additional BAO information by 
exploiting higher-order N-point statistics. 

Modeling reconstructions in LPT and expanding all fields consistently up to 0(5%) shows that Eulerian and La¬ 
grangian reconstructions change the density in exactly the same way (if a new correction term is included that arises 
by modeling the displacement field of clustered catalogs such that it is evaluated at Eulerian rather than Lagrangian 
locations, which matters beyond linear order). We confirmed in simulations that density changes due to Eulerian 
and Lagrangian reconstructions look nearly the same in 2D density slices, see Fig. [lO] Given that Eulerian and 
Lagrangian algorithms agree at the field level in theory and simulations and that they enhance BAO information at 
the power spectrum level in a very similar way, they can be regarded as Eulerian and Lagrangian incarnations of the 
same underlying reconstruction principles. By making this connection one can argue that the standard Lagrangian 
LGS reconstruction is partially so successful because it automatically combines BAO information from specific well- 
motivated 2-, 3- and 4-point statistics of the unreconstructed density. This provides a novel theoretical argument for 
the robustness and success of the standard LGS reconstruction algorithm. 

Using specific splits of the mass density in simulations, we show that most of the additional BAO signal-to-noise from 
Eulerian reconstructions is due to 3-point statistics, while 4-point statistics add very little. The 3-point contribution 
can be split further into parts of type (<5o |<5q) and (<5 q|< 5§), where <5o is the linear density. Our simulations show that 
the 13 part of the 3-point contribution is responsible for sharpening the BAO wiggles, whereas the 22 part leads to 
shifts and damping of the wiggles, in agreement with previous literature HE- 

The specific 3-point statistics that enter Eulerian reconstructions come in the form of cross-spectra of the unrecon- 
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structed density with the squared unreconstructed density, 8 2 , a shift term, s ■ VS, and a tidal term, K 2 . The same 
cross-spectra also arise when considering maximum-likelihood estimators for the amplitudes of the components of the 
tree-level DM bispectrum [Mj. This suggests that the bispectrum of the unreconstructed density does not contain 
much additional BAO information beyond that already recovered by reconstruction algorithms (i.e. the additional 
BAO information that the bispectrum could in principle yield is likely very correlated with the additional BAO in¬ 
formation in the reconstructed power spectrum compared to the unreconstructed power spectrum). This intuitive 
expectation should be tested more rigorously and quantitatively in the future. 

An important limitation of the results presented here is that we restricted ourselves to dark matter in real space. 
We leave the important extension to galaxies in redslrift space for future work. Our findings can also be extended 
in many other directions. For example, it would be interesting to study configuration space correlation functions in 
simulations and theory, and to include higher order modeling corrections. Ideas from previous studies mmm that 
aim at optimizing the standard reconstruction algorithm, e.g. by improved weightings or iterative reconstructions, 
could also be applied to each of our new algorithms. Our results could also be useful for improving combinations of 
BAO analyses (fitting only for the BAO scale with templates) with redshift space distortion analyses (using the full 
shape of the redshift-space 2-point function), because the covariance between these analyses could be expressed in 
terms of covariances between 2-, 3- and 4-point functions that can be modeled perturbatively. Ultimately it would be 
exciting to apply the new algorithms to real data. 


ACKNOWLEDGEMENTS 

We are grateful to Uros Seljak for initial collaboration and many helpful discussions. We also thank Nikhil Pad- 
manabhan and Martin White for useful discussions, and Martin White for providing the RunPB TreePM N-body 
simulations that were used for parts of this paper. This research used resources of the National Energy Research 
Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of 
the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. 


Appendix A: Lagrangian Modeling 

This appendix models Lagrangian reconstruction algorithms using the Lagrangian perspective of structure formation 
following [15] ■ We start by modeling LGS reconstruction non-perturbatively in terms of the displacement field. To 
simplify this, we then review Lagrangian Perturbation Theory and use it to model LGS reconstruction perturbatively 
up to second order in the linear density. Finally, we model all other Lagrangian reconstruction algorithms mentioned 
in the main text. Throughout this appendix we include a second order correction term that arises when assuming 
that clustered catalogs are displaced from Eulerian rather than Lagrangian positions. 


1. Nonperturbative Lagrangian Modeling of Lagrangian Growth-Shift Reconstruction 


In the Lagrangian picture of structure formation, initial positions q in Lagrangian space are displaced by \F(q) to 
obtain Eulerian positions x, 


x = q+ 'F(q). 


This can be related to the mass overdensity by m 


<5 th (x)= / d d qMx-q-*(q))-l, 


(Al) 


(A2) 


where <5p is the Dirac delta and we use superscript ‘th’ to denote theoretical mass densities. The Fourier transform is 


(5 th (k) = J d 3 x e" lk x J th (x) = J d 3 q e _ik q (e~ ik * (q) - l) . 


(A3) 


The LGS-reconstructed density is obtained as follows. First, the negative Zeldovich displacement s is calculated 
from the non-linear density with Eq. (12). This is used to displace the original density, yielding the ‘displaced’ field 
Sd which can be modeled by 


5‘ h (k) = J d 3 q e" lk ' q (g-A-[*(q)+s(x)] _ ^ . 


(A4) 
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The theory displacement is evaluated at the initial Lagrangian position q, whereas the negative Zeldovich dis¬ 
placement s obtained from the observations is evaluated at the final Eulerian location x. Instead of s(x) previous 
literature uses s(q) which misses a correction at second order in the fields (this correction is investigated further with 
simulations at the end of Appendix Jb|. Next, a uniform or random sample of particles is displaced by the same vector 
s evaluated at the uniform positions q to get the ‘shifted’ field S s , which can be modeled by 

<5‘ h (k) = J d 3 q e _ik q (V ,;ks(q) - l) . (A5) 

The difference of these two fields gives the reconstructed density for the standard LGS reconstruction algorithm, 

= (A6) 


which can be modeled by 


CGS h ( k ) = J d3 q e_ik q (V ik ' [ * (q)+s(x)1 - e“ ik ' s(q) ) . 


This can be rearranged nonperturbatively t< 


rrec.th _ rth , erec,(a) rrec,(b) 

°LGS — 0 T °LGS T °LGS > 


where 5 th models the nonlinear density before reconstruction and we defined 

<S[2s a) ( k ) = J d 3 q e" lk q ^ e -ik-[*(q)+8(q)l _ g-ik .^(q) _ g-ik-s(q) + ^ _ 


(A7) 


(A9) 


(A10) 


An additional correction term (b) arises because the clustered catalog displacement is evaluated at Eulerian positions, 
s(x), instead of Lagrangian positions, s(q), 


<^LGS b) ( k ) = J d 3 qe- k q G lM(q ) ( 


e -ik-s(x) _ e -ik-s(q) 


Expressing the exponentials in Eq. (A10) as power series gives 

„ oo n— 1 


c-rec,(a) 

5 lgs 


(k) = / d“q e-"“>££ 


(— i) n ( n 


7.! 


[k.*(q)r[k.s(q)] r 

'll.I \ III. / 

n =2 m—1 

where only terms mixing SI/ and s survive (m > 1 and n — m> 1). Similarly, if we write 

s(x) = s(q) + H(q) 


(All) 


(A12) 


for some function H (discussed in more detail later, see Eq. (A28)), then Eq. (All) becomes 


rrec,(b) 

3 lgs 


(k) = / d 3 q e -'k q e -*k *(q) Y 


HY 


n =1 m—1 


[k • H(q)] m [k • s(q)] T 


(A13) 


All expressions above are fully nonperturbative in the displacement field. They can be simplified further in La¬ 
grangian Perturbation Theory (LPT) which we briefly review in the next section. 


2. Review of Lagrangian Perturbation Theory (LPT) 

In LPT the displacement SI/ is expanded perturbatively in powers of the linear density contrast So, 

= SI/l 1 ) + ^ (2) -I-, (A14) 


12 This can be seen by adding and subtracting 5 th given by Eq. (A31, which yields 

-ik-q 


c-rec.in/i \ r-t 

^LGS ( k ) =<5 


'(k) + j d 3 q, 


+ 1 


(AS) 
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where 


= h / n 


i= 1 


d 3 k, 

_(2tt) 3 


(2^) 3 fe(k 1 + • • • + k„ - k) LW(k!,..., k n ; k) 5 0 (k!) • • • <5 0 (k„). 


L ( 1 ) (k) = 


The first and second order kernels are (see e.g. [291 [431446] ) 

k 

¥' 

L ( 2 ) (ki,k 2 ;k) = ^ 1 - P 2 (ki • k; 

where k = k/fc and P 2 is the l = 2 Legendre polynomial. 

The density contrast (A3) corresponding to this perturbative displacement field is 


^ th (k) = ^LP T (k) + <^Lp T (k) H 


where ^lpt = and 


A 2 ) 

J LPT 


(k) = / d 3 qe- ikq 


-ik- <ld 2 )(q) - i (k- ^ (1) (q))' 


(A15) 

(A16) 

(A17) 

(A18) 

(A19) 


noting that ^"^(q) are in Lagrangian configuration space. To simplify the comparison with Eulerian SPT, we write 
the second order held in terms of the kernel .Ep PT defined by 


#r(k)= [ i' , 2 LPT (k 1 ,k 2 )5 0 (ki)<So(k 2 ). 
J k, 


(A20) 

where we used the shorthand integral notation of Eq. |3j). The FJ^^ kernel has components from VJ/PlvI/l 1 ) in 

Eq. (A19) and i^ PT from i.e. 


LPT 


witl £3 


| 77. 

*2 ~ * 2,11 + ^ 2,2 


E 2 LPT ( kl) k 2 )-? + i (| + |) kl -k 2 
f ppt ( k 1 > k a ) = i-ip 2 (k 1 .k 3 ). 


gPaCkr-ka), 


This agrees with the F 2 kernel from Eulerian SPT [! 

F 2 LPT (k 1 ; k 2 ) = F 2 SPT (ki,k 2 ) = 


17 1 (k x 
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fc 2 k\ 


k 1 -k 2 + -P 2 (k 1 -k 2 ). 


(A21) 

(A23) 

(A24) 

(A25) 


The E^ PT contribution from L^ 2 ) does not change the shift term, but it changes the amplitude of the growth term 
by 21% and that of the tidal contribution by 43%. 


3. LGS Reconstruction in 2LPT 


Working in LPT and keeping terms up to second order in the linear density S 0 , the contribution (a) of Eq. (A12) 
to the reconstructed density becomes 

CGs a) ( k ) = - [ [k ■ L ( 1 ) (ki)l [k ■ L«(k 2 )l M / A(fc 2 )<5o(k 1 )<5 0 (k 2 ), (A26) 


/ k 7 L 


13 Note that for k = ki + k 2 we have 

[k- L (1) (ki)] [k-L«(k 2 )] = | + ki k 2 + ^P 2 (ki • k 2 ). (A22) 

Also note that since -^^(ki, —ki) = 0 and 1 ,—ki) = 0 the corresponding contributions to ^lpt ^ ave zero spatial average, 

(^lpt^W) = 0 and <4 pt( x )) = °- 
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where we wrote Sl'(q) and s(q) in Fourier space in terms of L/ 1 ). From Eqs. (A22) and (A23) we get m 

CGs a) ( k ) = ~J 2F^(k 1 ,k 2 )W R (k 2 )6 0 (k 1 )8 Q (k 2 ). (A27) 


Next, we calculate at second order in LPT the correction term (b) given by Eq. ( |A13| ) (coming from evaluating 
clustered catalog displacements at Eulerian rather than Lagrangian locations, see also Appendix [B]). At second order 
we have 


s(x) « s(q) + ('E'(q) • V)s(q), 


(A28) 


so that Eq. (A13) becomes at second order 

€cs b) (k) ~~ij d 3 q e- ik q e- ik s ( q )k • (*(q) • V)s(q). 


Then 


Ccs b) (k) = ^ (k • L«(k 2 )) (k 2 ■ LW(ki)) W R (k 2 )8 0 (ki)d 0 (k 2 ) 

w^(fc2)<y 0 (ki)Jo(k 2 ). 


/k i 


X + T“ki • k 2 + ^P 2 (ki • k 2 ) 

o o 


(A29) 

(A30) 

(A31) 


Adding Eq. (A31) to (A27) cancels the tidal part and one of the two shift terms of . so that the model for the 

total reconstructed density becomes simply 


Ccs th (k) = <5 th (k)- 


f* 

jfe, - 

L 

1 + tt ki • k 2 

L fc 2 J 


bbA(fc 2 )<5o(ki)<5o(k2), 


(A32) 


which agrees with Eq. (45) in the main text and corresponds to Eq. (46) in configuration space. Thus, in this LPT 


picture, standard LGS reconstruction subtracts growth and shift terms from the nonlinear density <5. 

The second-order contribution to the density change caused by reconstruction came from terms with kernels of 
the type (L^) 2 but no L^ kernel contributed, which means that this reconstruction cannot fully reverse second- 
order non-linearities ns. We can generalize this statement using Eq. ( |A12[ ): At n-th order LPT, part (a) of the 
reconstruction does not have an explicit contribution from the kernel , but only products of lower-order kernels 
L^i)... L( r <>) (with s > 2 and JA r* = n). For example, at third order only kernels of the type (L*A) 3 and (iP^) 2 lS 2 ' > 
contribute but not L^ 3 \ The same statement is true for the (b) part of the reconstruction given by Eq. ( |A13 Ip*] 
Consequently, the standard LGS reconstruction cannot fully remove any of the ‘pure’ L^ contributions to the 
nonlinear density, unless products of lower-order kernels conspire to agree with L^ n \ It can of course still improve 
BAO information by removing ‘non-pure’ non-linearities, e.g. of type (L^) 2 , and this turns out to be very successful 
in practice. 


4. 2LPT Models for 6 Lagrangian Reconstructions 


We introduced 6 possible Lagrangian reconstruction algorithms in Section m At first order, they all agree with 
the linear density, but at higher order they differ from each other. This appendix models these algorithms in 2LPT. 

The four relevant densities for the reconstructions are clustered catalogs displaced by s or —s, denoted by 5,i [s] and 
8d\— s ], an d ran dom catalogs shifted by s or —s, denoted by <5 s [s] and <5 S [—s]. We model them in 2LPT by expanding 
Eqs. (A4) and (A5) consistently at 0(6$). For clustered catalogs, the second order part gives 


4 2, [±s](k) = 


{[It w R (k)\ F 2 ( k 1; k 2 ) T w R (k 2 ) 


1 T y-^-ki • k 2 

k 2 


W R (fci) W R ( k 2 ) Ef if (ki, k 2 


d 0 (k 1 )<5 0 (k 2 ). 

(A33) 


14 The only possibility where a summand in Eq. jA13| is not a product of at least two fields is if m = n = 1, in which case there is a 
contribution going like H^ n ). However, Taylor expanding H(q) = s(q + — s(q) in ^ implies that this is always a product of at least 

two fields that depend on L kernels (schematically, suppressing derivatives and vectors, H ~ 'I's + 2 s + + •••). 
























25 


Here we used s(x) = s(q) + H(q), where lT 2 )(q) = (^(^(q) ■ V)sA)(q) and —ik • H^ 2 ^(k) is given by Eq. (A31). 
Also, note that the negative Zeldovich displacement s(q) computed from the clustered catalog is 


s (2) (k) = -iLW(k)W R (k) J F 2 ( k 1 ,k 2 )5 0 (k 1 )5 0 (k 2 ) 

at second order. For the shifted randoms, the second order part is 

^ 2) [±s](k) = f { TW / fi(fc)F 2 (k 1 ,k 2 ) + fE fl (fc 1 )fE«(fc 2 )F 2 L P 1 T (k 1 ,k 2 )}«5 0 (k 1 )^ 0 (k 2 ). 


The reconstructed densities of Eqs. (|38|)-(|43|) up to second order are then 

l + M} <5o( k 1 ) <5o(k2 ), 


S d [s] - 6 s [s] = <5 0 (k) + 

f {T 2 (k 1; k 2 ) 

{<5d[s] + Ad [—s]} = <5o(k) + 

f {^(k 1; k 2 ) 

Ad [s] + A s [—s] = <5o(k) + 

/ {T 2 (k 1; k 2 ) 

5 s [s] + <5<j[—s] = ^o(k) + 

f {*2(ki,k 2 ) 

{As [s] + A s [—s]} = <5o(k) + 

/»* 

1 {E 2 (ki,k 2 ) 

Ad[-s] -A s [—s] =(5 0 (k) + ^ 

1 {T 2 (k 1; k 2 ) 


(A34) 


(A35) 


(A36) 




|^o(ki)<5o(k 2 ), 


(A37) 




|^o(ki)<5o(k2), 


where p, = ki k 2 . From the six possibilities, only combinations (A36l and (A37l suppress nonlinear growth and shift. 


Appendix B: Slice Comparisons 

To compare reconstruction algorithms at the field level and as an illustrative sanity check of their 2LPT models, 
we write 


C>rec = 6 - (5 - 5 iec ). 


(Bl) 


Here, S — £) rec is the nonlinear excess density that is removed from the full nonlinear density when applying recon¬ 
struction. Fig. [To] shows 2D slices of this excess density for the three Lagrangian reconstruction algorithms (top) 
in comparison with the three corresponding Eulerian algorithms (bottom). For each column of the plot, upper and 
lower panels are rather similar, i.e Eulerian and Lagrangian reconstructions change the density in a very similar way. 
This shows qualitatively that Eulerian and Lagrangian reconstructions agree with each other at the field level on the 
large scales that are highlighted by these slice plots. The agreement between upper and lower panels in Fig. [To] also 
illustrates that second order LPT provides an accurate description of the reconstruction algorithms on large scales. 
The different Lagrangian reconstructions shown in the different columns are also rather correlated, but they differ 
somewhat, which is expected since the 2LPT models of these methods differ from each other. 

Fig. [lT] shows some additional slice plots for reference and comparison. The density before reconstruction is a few 
times larger than typical changes due to reconstruction. The growth term 5r5 peaks at density peaks and troughs 
which are amplified by squaring the field. The magnitude of the shift term s • V<5 is maximal at locations where the 
large-scale Zeldovich displacement is aligned with the density gradient, which is most pronounced near density peaks. 
Since growth and shift terms have opposite signs near density peaks, their combination cancels the largest peaks to 
leave a smoother combined field (see Fig. 10). 


Finally, we check if displaced clustered catalogs should be modeled by evaluating the displaceme nt fi eld s at Eulerian 
positions x rather than Lagrangian positions q, which leads to the correction term (b) in Eq. (A91. The panel on 
the right in Fig. |TT] shows the expected excess density for LGS reconstruction if that correction term is ignored and 
only term (a) is used in Eq. (A91. The same plot including the correction term (b) is shown in the bottom left panel 
of Fig. [To] These two panels should be compared with the top left panel of Fig. [lO] which shows the true excess 
density for the full LGS reconstruction. The impact of term (b) is not large but it does significantly smoothen some 
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FIG. 10. Comparison of Lagrangian (top) and Eulerian (bottom) reconstructions at the field level for growth-shift (left), random- 
random (middle) and F2 (right) algorithms. The panels show 2D slices of excess densities 8 — 8 lec that are removed from the 
density when performing reconstruction. The excess densities in the bottom panels are quadratic in the nonlinear density 
measured from simulations. They coincide with the excess densities expected from modeling the Lagrangian reconstructions 
in the upper panels with 2LPT, if the linear densities in Eqs. (461, (481 and (531 are replaced by the nonlinear ones. For LRR 
and ERR, twice the excess density is shown to enhance color contrast. We use R = 15h _1 Mpc Gaussian smoothing for the 
reconstructions. To highlight large scales, each final field is additionally smoothed externally with a R = 15h -1 Mpc Gaussian. 
The clustered catalog is obtained from a 1% random subsample of one RunPB N-body realization; the random catalog from 
equally many randomly distributed particles. 
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FIG. 11. Left panel: Nonlinear density before reconstruction (divided by 3 for better visibility). Middle panels: Individual 
growth and shift contributions to the Eulerian growth-shift reconstruction. Right panel: Expected excess density for growth- 
shift reconstruction if the correction term (b) in Eq. (A91 is ignored. The plot shows term (a) given by Eq. (A271 (multiplied 
by -1). All densities were obtained from the same 1% subsample as in Fig. 10 and are plotted in the same way. 
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of the peaks and troughs of the density (slightly difficult to see in the slices because colors are truncated) and leads 
to slightly better agreement with the excess density for the full LGS scheme. The slice plots thus suggest that the 
correction (b) should be included in the 2LPT modeling, although its effect is not large. 

It should also be noted that the Eulerian growth-shift reconstruction derived at the beginning of the paper from the 
continuity equation only agrees with the standard Lagrangian reconstruction in 2LPT if the correction (b) is included 
in the modeling of the latter. An Eulerian reconstruction scheme that would only resemble term (a) of Eq. (A91 yields 
a cumulative BAO signal-to-noise-squared of ~ 150 for fc max = 0.4/i/Mpc in our simulations. This is significantly less 
than that for the EGS method which includes the correction term (b) and resembles the BAO information of LGS 
reconstruction more closely. We take this as additional evidence that term (b) should be included. 


Appendix C: Newton-Raphson Method 


We introduced Eulerian EF2 reconstruction in Section m by subtracting from the nonlinear density the second 
order F% part, which was approximated to be quadratic in the observed rather than linear density. Here we take a more 
formal approach to motivate this algorithm using the Newton-Raphson iteration to find a linear density that generates 
a given observed nonlinear density (at second order). Note that our approach differs from the (more complicated) 
iteration scheme proposed in [T 8 ] , which is based on maximizing the likelihood of the displacement field under certain 
assumptions. 

In Eulerian Standard Perturbation Theory (SPT), the theoretical nonlinear DM density is modeled in terms 
of the linear density Sq as 

4l(x) = <So(x) + F 2 [(So](x) H - , 

where the second-order perturbation is 


(Cl) 


^N(x) = ^o(x) - ^o( x ) • V£o(x) + ^y-ffo(x), 


(C2) 


and the linear displacement tko and linear tidal term are defined in terms of So by Eqs. (24)-(26). We neglect higher- 


order terms like F 3 [<5q]( x) in Eq. ( |C1| ). Given a fixed observed non-linear DM density <5 0 bs, our goal is to estimate the 
linear density So that generated this observed density, i.e. we want to find So such that 

f[S 0 ] = ^0 + A 2 [(>o] — £>obs = 0, (C3) 

i.e. we want to find roots of the functional / with respect to So- This can be achieved with the Newton-Raphson 
iteration method, where the (N + l)-th iteration step is 

r +1) (k) = 4 iV) (k) - J M N) W, (C4) 


where superscripts label iteration steps throughout this section. /' in Eq. (C4) is the functional derivative of /. This 
is a linear operator that can be written in matrix form as 

d/[*o(k)] 




Ojk.k' — 


= Mk - k') + 2 E 2 (k', k - k')<S 0 (k - k'). 


d<5 0 (k') 

For S 0 < 1, the inverse of the derivative is then 

(f’[S oDki = Mk - k ') - 2 F 2 (k', k - k%(k - k') + O(S 2 0 ). 
The Newton-Raphson iteration scheme thus becomes 

d 3 k' 


(C5) 


(C6) 


<T +1) (k) = 5 obs (k) - EM^Kk) + 2 


(2 7 T ) 3 


F 2 (k',k-k , )4 Af) (k-k') ^(kV.Wk') TO)^) 3 )- (C7) 


If we start the iteration with 5^ = 0, then the first iteration step gives 5^ = 5 0 bs) he. the linear density is 


approximated by the non-linear observed density! 33 ] The second iteration step gives 


4 2) ( x )=^obs-F2[<5obs]+0(5 3 bs ) 

= ^obs(x) - g^(x) + ^obs(x) ■ V<5obs(x) - ^-A'obsW + ^(^obs) 


(C8) 

(C9) 


15 


We could as well start the iteration directly with = 5 0 b s ; then the N-th step would agree with the (N+l)-th step of the iteration 
started with 5^ = 0. 
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^max = 0.4/l/MpC 



FIG. 12. Cumulative BAO signal-to-noise for varying Gaussian smoothing scale R, using scales up to fc max = 0.15/i/Mpc (left) 
or fc max = 0.4/i/Mpc (right). Note that the shot noise in our setup is unrealistically small, h _1 = 31Mpc 3 //i 3 . 


This agrees with the EF2 algorithm from Section III (up to smoothing, which we ignored here for simplicity). 


Appendix D: Numerical Analysis Details 
1. Simulation Convergence Tests 

To check convergence of the FastPM simulations for the quantities and plots relevant to this paper, we perform 
several convergence tests. First, we generate a single FastPM realization that matches the cosmology and phases of 
the initial conditions of one of the TreePM RunPB simulations. We then measure all auto- and cross-spectra between 
<5, A 2 , HI • VS and K 2 , corresponding to 2-, 3- and 4-point statistics required for Eulerian reconstructions. All spectra 
from FastPM are slightly lower than the corresponding TreePM RunPB spectra over the whole k range, and the 
deviations are smallest at low k. The density power spectrum (2-point) is less than 0.8% low over the whole k range, 
3-point spectra like (<5 2 |d) differ by less than 1 % at k < 0.2/i/Mpc and by at most 1.5% at higher k. 4-point spectra 
like (<5 2 1<5 2 ) differ by less than 1.5% at k < 0.2/i/Mpc and by at most 2% at higher k. This is true for full DM samples 
as well as 1% subsamples (with the same particles selected from the RunPB and FastPM simulation). In summary, 
FastPM simulations and TreePM RunPB simulations agree at the 2%-level. 

Since BAO wiggles are only percent-level fluctuations on top of the broadband shape of the power spectrum, even 
differences at the 2% level could potentially be problematic for BAO studies. We expect however that broadband 
systematic inaccuracies in the simulations should cancel out when forming differences between wiggle and nowiggle 
simulations, as long as these systematics are present in wiggle and nowiggle simulations, which is a reasonably 
assumption. To test this more quantitatively we vary the accuracy settings of the FastPM code and check if any 
results depend on them: First, we run the same simulations with only 20 time steps and starting redshift z,- t = 9 
(instead of 80 time steps and Zi = 99 used for the fiducial simulations). Second, we also run them with 80 time 
steps and Zi = 9. We generated every plot of this paper based on FastPM simulations for the three different accuracy 
settings. All plots agree almost perfectly between the three settings, and changes are so small that they cannot 
be seen by eye. Therefore, our results do not depend on the accuracy settings of the FastPM code. This provides 
additional evidence that the simulations are converged for the purposes of this paper. For all plots that require 
nowiggle simulations, we will only show plots for the fiducial FastPM simulations with 80 time steps and Zi = 99 
because they should be most accurate. For plots that do not require nowiggle simulations we use RunPB TreePM 
simulations instead. 


2. Dependence on Reconstruction Smoothing Scale 

All numerical results in the main text assumed a fiducial Gaussian smoothing scale of R = 15Mpc/ft. for the 
reconstruction algorithms, which was found to be optimal for the analysis in [13 Q3J. Here we briefly discuss how our 
results depend on this choice. Fig. |l2| shows the cumulative BAO signal-to-noise for different smoothing scales R. For 
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fc max = 0.15/i/Mpc, Eulerian and Lagrangian reconstructions perform very similarly independent of smoothing scale, 
which is expected because the 2LPT models of the algorithms agree. If smaller scales are included, k m ax = 0.4Mpc//i, 
Eulerian and Lagrangian reconstructions still perform similarly for large smoothing scale, R > 15Mpc /h. In particular, 
EGS, EF2 and LGS all have the same performance for R = 30Mpc//i. Towards smaller smoothing scales, especially for 
R ~ 5Mpc//i, the LGS and LF2 algorithms improve significantly, whereas their Eulerian counterparts EGS and EF2 
do not improve as much. This shows the limitations of the Eulerian methods on small scales k > 0.15/i/Mpc and small 
smoothing scales R ~ 5Mpc//i, where Lagrangian algorithms seem to work somewhat better; EGS misses up to 13% of 
the total BAO signal-to-noise of LGS in the most extreme case (for R = 2.5Mpc /h and fc max = 0.4/i/Mpc). However, 
it is important to note that the shot noise of our 1% DM subsamples is unrealistically small, n -1 = 31Mpc 1 2 3 4 5 6 7 8 9 //i 3 . For 
more realistic, higher shot noise levels, the shot noise power starts to dominate on larger scales so that reconstruction 
performance is expected to be optimal for R ~ 15Mpc/!i and decrease for smaller smoothing scales [TO]. We also 
checked that our results do not qualitatively change if the Gaussian smoothing kernel is replaced by a steeper logistic 
smoothing function Wn-s{k) = l/{l+exp[S'(fc — 1/R)}}, where we tried R = 5Mpc /h and R = 10Mpc//i with S = 100. 


3. Gaussian Covariance 


Throughout this paper we assumed diagonal Gaussian covariances given by Eq. (601. We briefly argue here why 


this is a reasonable approximation for the purposes of estimating the BAO signal-to-noise. 

The non-Gaussian off-diagonal corrections to the power spectrum covariance that become important in the nonlinear 
regime can be modeled by m 


co v(P(fc),P(A/)) = (P(k))(P(k')) 


&k,k 


N, 


modes 


(k) 


CNG 


(Dl) 


where the non-Gaussian correction cng is independent of k and k' . It arises from the non-Gaussian coupling of large 
and small scale modes that leads to a broadband up- or down-shift of the small-scale power spectrum depending on 
the specific realization of large-scale modes in any given realization (see e.g. [34)). In other words, in some realizations 
the ratio of wiggle over nowiggle power spectra will be high over the whole k range, while in other realizations 
that ratio will be low over the whole fc-range if large-scale modes happen to be different (we confirmed this effect 
in our simulations). In any given realization, the broadband rescaling of the power spectrum will be absorbed by 
the broadband part when fitting for BAO wiggles on top of some generic broadband function (as done in practice). 
Therefore the non-Gaussian correction cng in Eq- (Dl) does not affect the BAO signal-to-noise and can be ignored 
for the purposes of this paper. 

We confirmed in simulations that the variance of spectra estimated from the mode-to-mode scatter in any single 
realization is in excellent agreement with the Gaussian expectation. The variance estimated from the scatter between 
our FastPM realizations 1 and 2 also agrees with this, but it is significantly larger when estimated from the scatter 
between realizations 1 and 3 or 2 and 3. The reason is that our realization 3 has a higher overall amplitude over the 
whole k range than the other two realizations, which is compatible with the arguments above and does not affect the 
significance of BAO wiggles in any given realization. 
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